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Utilizing large-scale Monte-Carlo simulations, we investigate an unconventional two- component 
classical plasma in two dimensions which controls the behavior of the norms and overlaps of the 
quantum-mechanical wavefunctions of Ising-type quantum Hall states. The plasma differs funda- 
mentally from that which is associated with the two-dimensional XY model and Abelian fractional 
quantum Hall states. We find that this unconventional plasma undergoes a Berezinskii-Kosterlitz- 
Thouless phase transition from an insulator to a metal. The parameter values corresponding to 
Ising-type quantum Hall states lie on the metallic side of this transition. This result verifies the 
required properties of the unconventional plasma used to demonstrate that Ising-type quantum Hall 
states possess quasiparticles with no n- Abelian braiding statistics. 
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I. INTRODUCTION 

Key properties of physical systems can sometimes be 
understood by mapping them to seemingly unrelated 
ones. A powerful example of this was provided by Laugh- 
lin, who observed that the squared norm of his v — 1/M 
fractional quantum Hall trial wavefunction 
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(where Zi = Xi + iyi is a complex coordinate in the two- 
dimensional plane) could be expressed as the Boltzmann 
weight of a two-dimensional one-component plasmsW 
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where 

vite) = -Qi f> \zi - zj\ + H J2 M 2 (3) 
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and Qi/T = 2M. This mapping allows properties such 
as quasiparticle charge and braiding statistics to be de- 
termined by appealing to the known properties of a one- 
component plasma. 

Recently, a similar plasma mapping was established 2 
for Ising-type quantum Hall states, such as the Moore- 
Read (MR) 3 , anti-Pfaffian 4 5 , and Bonderson- Slingerland 
(BS) hierarchy 6 states, which are likely candidates to de- 
scribe Hall plateaus in the second Landau level, in par- 
ticular at filling fraction v = 5/2P"^. In this case, the 
mapping is to a two-dimensional (2D) two-component 
plasma, where the two species of particles, w and z, 



carry not only different values of charge, but also interact 
through two different interactions, both of the Coulomb 
form, so the potential energy is: 

V{z i] w a ) = V 1 {z i ) + V 2 {z i] w a ), (4) 

N N 
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+ Q 2 ^ln|^-w a |, (5) 
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where Q 2 /T = 3. The z-particles interact with each 
other through the first Coulomb-like interaction, Vi(zi), 
given in Eq. |3| (and so does not depend on the w a 
coordinates). Moreover, the z-particles interact with 
each other and with the ^-particles through the second 
Coulomb-like interaction, through which the ^-particles 
also interact with each other, according to ^2(^5^)5 
given in Eq. Note that Vi(zi \ w a ) is the 2D Coulomb 
potential of the usual two-component plasma (where the 
two species carry charge Q2 and — Q2, respectively). 

The z-particles carry charge Qi for the first interaction 
and charge Q2 for the second interaction. The u>-particles 
carry charge for the first interaction and charge —Q2 for 
the second interaction. For a plasma with N particles of 
each species, neutrality is satisfied using a uniform back- 
ground density of type 1 charge, as in the second term in 
Eq. (|3|. This unconventional plasma may be considered 
as an ordinary neutral two-component gas with positive 
and negative charges of magnitude Q 2l where the posi- 
tive charges are given an additional charge of Qi that is 
only felt by the other positive charges and not the nega- 
tive charges. An illustration of the interactions between 
the two species in the system is shown in Fig. [T] 
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Figure 1: (Color online) Illustration of interactions between 
the particles in the 2D system. The if -particles only in- 
teract by the second Coulomb-like interaction with charge 
— Q2, whereas the z-particles carry charge Qi for the first 
Coulomb-like interaction and Qi for the second Coulomb-like 
interaction. Thus, the intraspecies interaction among the w- 
particles, shown in (a), and the interspecies interaction be- 
tween w and z-particles, shown in (b), are given by Qi only, 
whereas the intraspecies interaction among the z-particles, 
shown in (c), are determined by Q\ in addition to Qi. Interac- 
tions between the z-particles and the neutralizing background 
are omitted from the figure. 



We are thus led to consider a class of unconventional 
plasmas parametrized by Q\/T and Q\jT . As mentioned 
above, for MR Ising-type states with filling v = 1/M, 
the relevant values are Q\/T = 2M and Q\jT = 3. In 
this plasma mapping, the Z{ particles in the plasma cor- 
respond to the electrons in the MR wavefunctions and 
the w a particles correspond to screening operators (ficti- 
tious particles). The case Qi = 0, Q\jT = 3 is relevant 
for the plasma mapping^ of 2D chiral p-wave supercon- 
ductors^. We note that whenever Qi =0, our model 
is a special case of the well-known 2D two-component 
plasma of equal and opposite charges^HHI. The screen- 
ing properties of multi-component 2D plasmas with mul- 
tiple Coulomb interactions of this kind are also impor- 
tant for other physical systems, such as rotating multi- 
component Bose-Einstein condensates with interspecies 
current-current (Andreev-Bashkin) interactio n 16 * 17 * and 
some mult i- component superconducting systems^MD i n 
these systems the screening properties and phase transi- 
tions determine superfluid and rotational responses. 

In this paper, we fix temperature to T = 1 and consider 
the two most significant values of Qi, namely Qi = 0, 2. 
We investigate the screening and phase transition prop- 
erties of these plasmas as a function of varying Q2 by 
performing a large-scale Monte Carlo simulation. Here, a 
" screening phase" means that the system has a screening 
length which is finite, and exponentially decaying effec- 
tive interactions. A system with logarithmic effective in- 



teractions is one where screening is defined to be absent. 
As a first check, we reproduce the well-known result that, 
for Qi = 0, there is a Berezinskii-Kosterlitz-Thouless 
(BKT) phase transition at Q\ = Q\ c « 4, as expected 
for a 2D two-component plasma of equal and opposite 
charges. For Q\ < the charges are unbound and the 
plasma screens, but for Q\ > Q\^ the charges are bound 
into dipoles and the interaction is not screened. Thus, for 
Q\ = 3, the value relevant to 2D chiral p-wave supercon- 
ductors, the plasma screens. For Qi = 2, we again find a 
BKT phase transition at Q\ = Q\ c ~ 4, with a plasma 
screening phase for Q\ < Q\ Q . The first Coulomb-like 
interaction is deep within its screening phase and ap- 
pears to have a negligibly small effect on the screening of 
the second interaction. In both cases, the critical values 
Q\ c are obtained by a finite-size scaling fit of the Monte 
Carlo data to the BKT form. Our findings demonstrate 
that the unconventional plasma which occurs in the map- 
ping for both a chiral p-wave superconductor and the 
Ising-type quantum Hall states is clearly in the screening 
phase (for both types of Coulomb interaction) and hence 
allows one to discern the non-Abelian braiding properties 
of these states, as explained in Ref. [2j 

The outline of this paper is as follows. In the intro- 
ductory part of Section [TTJ we present the model for the 
unconventional plasma we will be studying in this paper. 
In Section [TlAj we connect this to the Ising-type of quan- 
tum Hall states. In Section II B[ we explain its connection 
to two-component, two dimensional, Bose-Einstein con- 
densates. In Section III A[ we present a formulation of 



the model on a sphere. In Section IIIB we give details 



of the Monte-Carlo simulations, and in Section III C[ we 
present our results for the screening properties, as well 
as our findings for the character of phase transition be- 
tween the dielectric non-screening phase and the metallic 
screening phase. In Section |IV[ we present our conclu- 
sions. Technical details on the derivation of a general- 
ized dielectric constant is given in Appendix [A] In Ap- 
pendix [Bj we give a derivation of a relevant higher order 
response function that we use to characterize the metal- 
insulator transition. In Appendix[Cj we present technical 
details on the finite-size scaling we have used. 



II. MODEL 



The canonical partition function of the unconventional 
plasma is written 



(6) 
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where the potential energy V is given by the 2D Coulom- 
bic interactions 

N 

V = Ql ^2 V ™w(\Wa-Wb\) 
a<b=l 

N 

+ (Qi + Ql) J2 vzz{\*i-*j\) 

i<j=l 

N 

+ V ™(\ Z i ~ W «D + ^,BG- (7) 
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Similar to the study of the 2D two-component neutral 
Coulomb gaPHmil we introduce a short-range hard- 
core repulsion between all charges in the system. Treat- 
ing all charges as hard disks with the same diameter d 
that limits the range of the hard-core repulsion, the in- 
teraction between charges of the same species is 

^(W) = ^(W) = {-ln|r|, \r\H ^ 

and the interaction between charges of different species 
is 

»~(M> = { In |r|, jrj ; d. W 

In Eq. Q, w a are position vectors for the particles of 
component w, and are position vectors for the particles 
of component z. To ensure neutrality, the term V Zi bg in- 
cludes the interaction of the Qi charges of type 1 for the 
z-particles with a neutralizing background charge den- 
sity. In Ref. [2j this background is a uniform negatively 
charged 2D disk with charge density qf G = — NQi/A, 
where N/A 1/2ttM, that yields 

1 N 

^bg—EN 2 ' ( 10 ) 

2=1 

The particle-background and the background- 
background interaction also yields uninteresting constant 
terms, that are disregarded in Eq. 

We note that when Qi = we have the 2D two- 
component neutral Coulomb plasma, wh ich is well- 
studied both analytically and numericall} ^ 15 1 25 ^ . At 
low dipole density, this system will undergo a BKT tran- 
sition, which is a charge-unbinding transition from a 
low-temperature state where charges of opposite signs 
form tightly bound dipoles to a high-temperature state 
in which a finite fraction of charges are not bound in 
dipoles, but rather form a metallic state. In the low- 
temperature phase, this Coulomb gas is an insulator and 
the dielectric constant e (see for instance Refs. I26|31|32l 
and Appendix [A| for a formal definition of e), is finite. In 
the high-temperature phase, the existence of free charges, 
yields a conductive gas with an infinite value of e. At 
the critical temperature T c , when tightly bound dipoles 



starts to unbind, there is a universal jump in the inverse 
dielectric constant from a non-zero value in the insulating 
phase to zero in the metallic phase, 

-i_ f 4T C , T->T C - ( . 

6 -\0, T->T+. (Uj 

The screening properties that follow are that the 
Coulomb gas is able to perfectly screen test charges in the 
metallic phase when there are free charges in the system, 
whereas there is no screening in the insulating dielectric 
phase. In this work, we will focus our attention on the 
low dipole density regime, so we will not go into details 
of the physics in the 2D two-component neutral Coulomb 
gas at higher densities. However, we note that when den- 
sity is increased, the critical point of t he BKT t ransition 
is shifted towards lower temperature d 14 l 15 l 28 l 29 i 

Another well-studied case is when Q2 = 0, for which 
the model reduces to the 2D one-component plasma 
(for the z-particles only). Early numerical studies of 
this system found a weak first-order melting transition 
at Qi/T ~ 140 from a state where the charges form 
a triangular lattice with quasi-long-range translational 
and long-range orientational order to a fluid plasma 
stat e 33 " 36 !. These results were, in a sense, contrasting 
with the defect-mediated melting theory of Kosterlitz- 
Thouless-Halperin-Nelson- Young (KTHNY) that pre- 
dicts melting from a solid to a liquid via two BKT- 
transitions and an intermediate hexatic phase with no 
tr anslat ional order and quasi-long-range orientational or- 
d e j 12 l 37tf 39[ Other studies of 2D melting point in favor 
of the KTHNY theorj^ES, suggesting that the nature 
of melting transition may depend on details in the in- 
teratomic potential, or that finite-size effects and lack 
of equilibration might lead to erroneous conclusions in 
earlier works. There are also studies that argue for 
the absence of a phase transition to a low-temperature 
solid phase in the 2D one-component plasma with repul- 
sive logarithmic interactions because the crystalline state 
would be unstable to proliferation of screened disclina- 
tions for any T > (pSHZl 



A. Ising-Type quantum Hall states 

The unconventional 2D two-component plasma studied 
here is mapped to inner products of trial wavefunctions 
for the MR quantum Hall states using conformal field 
theory (CFT) methods, as explained in Ref. |2j In par- 
ticular, this mapp ing utilizes the Coulomb gas descrip- 
tion of CFTsSSESI together with a procedure for replacing 
holomorphic-antiholomorphic pairs of contour integrals 
in screening charge operators for 2D integrals^^. 

The MR states' wavefunctions can be written as a 
product of correlation functions of fields from the Ising 
and U(l) CFTs. In particular, the MR ground-state 
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wavefunction for N electrons is 
9(z u ...,z N ) =Pf ' 1 



N 



i<j 
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where the PfafEan of an antisymmetric matrix A is given 

by 
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Here, 5V is the symmetric group, a is one of the permu- 
tation elements in Sn, and sgn(a) is the signature of a. 

The Pf (^ z ^ z . ^ portion of this wavefunction is produced 
from the correlation function of ip fields in the Ising CFT, 
while the Laughlin-type portion 
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is produced from the correlation function of vertex oper- 
ators in the U(l) CFT. 

The Laughlin-type portion of the MR wavefunctions 
can be mapped to charges of type 1, similar to Laughlin's 
plasma mapping. The mentioned CFT methods provide 
identities such as 
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(15) 



which allow the PfafEan portion of the MR wavefunc- 
tions to be mapped to charges of type 2. This allows 
one to write the norm of the MR ground-state wavefunc- 
tion as the partition function of the unconventional 2D 
two-component plasma of Eq. (|4|) 
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with Q\ = 2M and Q\ = 3. More generally, one can 
also construct a similar, but more complicated map- 
ping between inner products of wavefunctions of the MR 
states with quasiparticles, as explained in Ref. [2j In this 
case, the quasiparticles map to fixed "test" objects in 
the plasma that carry electric charge of type 1 and can 
carry both electric and magnetic charges of type 2 (and 
also changes the number of screening operators, i.e. w- 
particles in the plasma, to maintain neutrality). (The 
charges of type 1 and 2 carried by the quasiparticles are 



typically some fractions of the charges Qi and Q2 carried 
by the z-particles.) 



Strictly speaking, the right-hand-side of Eq. (15) is 
3 (since the integrand diverges as 



divergent for Q\ 



\w a — Zi\ as a w-particle approaches a ^-particle). It 
can be made well-defined (and equal to the left-hand- 
side) by replacing \w a — Zi\~ 3 with \w a — Zi\~ a , evalu- 
ating the integrals for a < 2 and then analytically- 
continuing a to 3. On the other hand, we regularize 



(13) the divergences of Eq. ( 16 ) in this paper by using a hard- 



core repulsion that forbids the particles from approaching 
each other closer than a distance d, i.e. replacing V in 
this expression with that of Eq. ([7]). It should not matter 
how we regularize the divergence in Eq. (16) as long as 



the probability for z-particles and u>-particles to sit right 
on top of each other has measure zero. As we will see 
in this paper, this is true for Q\ < Q\ c « 4, in which 
case the configurational entropy to be gained by having 
^-particles and ^-particles separate overcomes the energy 
gained by having them on top of each other. We refer to 
this as an "entropic barrier" for putting z-particles and 
u>-particles on top of each other. In contrast, in Eq. (15), 



where only the i^s are integrated over and the Z{ coor- 
dinates are fixed, regularization by a simple hard-core 
repulsion does not appear to be a suitable alternative 
to analytic continuation. In this case, since the Z{ coor- 
dinates are fixed, the entropic barrier is lower. Equiva- 
lently, there are fewer integrals to compensate the inverse 
powers. Thus, in Eq. (15), a simple hard-core cutoff will 



not reproduce the left-hand-side, and one must use the 
analytic continuation procedure described above. 



Two component rotating Bose-Einstein 
condensate in two dimensions 



In a rotating frame, a Bose-Einstein condensate in the 
London limit is described by the uniformly frustrated XY 
model, 



H 



d 2 r 



m 



V6>(r) - — 0(r) 



(17) 



where p = h 2 n/m for a condensate with mass m, phase 
9 and density n, and 0(r) = ft x r where ft = Qz is 
the angular velocity of the rotation. In 3D, this model 
is frequently used to describe the melting of vortex-line 
lattices in extreme type-II superconductors and super- 
fluid s 51 " 54 i By a duality transformation, the model in 
Eq. (flTl) can be rewritten in terms of vortex fields I to 
yielc 



H 



1 



d 2 q [Z(q) - (27T) 2 /5(q)] ± 
* ./ q 



[i(-q) - (27r) 2 /<H-q)] 



(18) 



where / = 2Vt/(j) is the vortex number density and 
0o = 27rh/m is the fundamental quantum unit of vor- 
ticity. This is a one-component 2D classical Coulomb 
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plasma where charges correspond to nonzero values in 
the vortex field l(r) and the quantity / now plays the 
role as the neutralizing background number density. 

Extending to two components, a model for a rotating 
two-component Bose-Einstein condensate with a generic 
Andreev-Bashkin drag interaction^^ reads 



H 



1 



. i=l,2 



hV6i 

mi'n t | fc) 

rrii 



^m 1 m 2 n d 



W(9i nvo 2 N 



mi 



m 2 



(19) 



where now m, n and # is given an index that denotes 
the component and is the drag density. This model 
has recently been studied in three dimension d 16 * 17 *. By 
a duality transformation, we arrive at the following 2D 
Coulomb plasma 



H=\J d 2 q [hid) - (27r) 2 /^(q)] 
x ^(-q) - (27r) 2 /^(-q)] , 



(20) 



where fa — 2fi/0o,i, </>o,i = 2ith/mi, U is the vortex field 
of component z, 



R = H 2 



mi 



Til 




1 



m 2 




and an implicit sum over repeated component indices 
z, j is assumed. By setting h = m^ = 1 such that 
/i = /2 = /? an d absorbing a factor 27T/3 in the den- 
sity coefficients, we see that the two-component Bose- 
Einstein condensate in Eq. (19) with n\ = 0, n 2 = Qi 
and = — Q\ corresponds to the unconventional two- 
component Coulomb plasma in Eq. ([7|. Thus, the un- 
conventional Coulomb plasma has a counterpart in a 
two-component Bose-Einstein condensate with a nega- 
tive non-dissipative drag interaction. However, note that 
in order to preserve a fixed number of charges when go- 
ing from the plasma description in Eq. ^ to the phase 
description in Eq. (19), we have to fix the number of 



vortices to only include rotationally induced vortices. In 
principle, in the BEC problem, the system can thermally 
excite vortex-antivortex pairs, but that process can be 
substantially suppressed by going beyond the phase only 
model in Eq. ( [l9| ) and introducing an additional energy 
penalty associated with vortex cores. 



III. MONTE-CARLO SIMULATIONS 
A. Considerations for a spherical surface 

Computer simulations of Coulomb interactions are 
generally difficult to perform due to the long-ranged na- 



ture of the interaction. Several techniques have been pre- 
sented to deal with the complications that arise 60 62 . We 
have performed large-scale Monte-Carlo simulations of 
the system described in Eqs. (|6| and (|7|) on a spherical 
surface. For other simulations on a spherical surface, see 
Refs. I14l34l43ft46l63l This may seem like a brute-force 
approach since the workload of the simulations scales as 
0(N 2 ). However, the benefit is that there are no bound- 
aries, the implementation is relatively easy, and there is 
no need to constrain the particles to move on a lattice. 
However, one must also be aware that simulation results 
may differ due to effects induced by topology. For in- 
stance, the triangular crystalline ground state of a 2D 
one-component plasma will necessarily include a num- 
ber of dislocations and disclinations on a sphere. These 
defects are not present in the ground state when the one- 
component plasma is located on the plane^El. 

We consider a sphere with radius R, with origin de- 
fined as the center of the sphere such that all particle 
position vectors w a and Zi are radial vectors with fixed 
magnitude R in three dimensions. The distance between 
the particles is measured along the chordlJJ 6 -^ 



■ r j ; | = 2i?sin 



where 



ipij = arccos(r^ • rj 



(22) 



(23) 



is the chord angle between the two particles at and rj 
with unit vectors r$ and fj, respectively. We may now 
rewrite the model in Eq. ^ on the surface of a unit 
sphere as 



V 



N 



N 



a<b=l 



N 



{Ql 



v zz (ii-ij) 

i<j=l 



(24) 



with interactions given by 



and 



-ln(l - Ti • ij), ipij > d/R, 



oo, ipij < d/R, 
ln(l -Ti-Tj), ipij > d/R. 



(25) 



(26) 



Note that the interaction V^bg m Eq. ([7]) between the 
neutralizing background and the excess charge of type 1 
becomes a constant term on the sphere, so we disregard 
it in Eq. ([M. 



The dimensionless density of particles on the sphere 
is given by the packing fraction r] = Ns/A where s = 
Asm 2 (d/4R) is the area of a hard disk of diameter d on 
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the sphere of area A = 4ttR 2 . In the simulation, we use 
a unit sphere with R = l. 

As explained in Appendix [A[ in order to account for 
screening properties when particles interact by two in- 
teractions simultaneously, we measure a general inverse 
dielectric constant, ej} x, given by 



(ai,a 2 ) 



C (ai,a 2 )' 



2aia 2 e 12 1 



l 2 t -22 ' 



where 



(27) 



(28) 



is a type specific inverse dielectric constant, a\ and a 2 
are type-dependent weights for the contributions of the 
different e~l (which are determined by the values of both 
types of charge carried by the test particles for which 
screening is being measured), and where Mi and M2 are 
the dipole moments for charges of type 1 and type 2, 
respectively, given by 



N 



Mi = Q 1 R^ f z i 



N 



N 



M 2 



(29) 



(30) 



Note that the type 2 inverse dielectric constant, e^ 1 , is 
the same dielectric constant as was used when studying 
the two- co mpon ent neutral Coulomb plasma on a spher- 
ical surface™^!. In addition to measuring the screening 
properties, the inverse dielectric constant may be used to 
identify the existence of a BKT-transition if it exhibits 
a universal discontinuous jump at the critical point, ac- 
cording to Eq. ( pT| . 

In addition to the inverse dielectri c con stant, we also 
measure the fourth-order modulus, " / 65 | 66 i This quan- 
tity may be used to verify a discontinuous jump in the 
inverse dielectric constant without making any a priori 
assumptions regarding the character of the phase transi- 
tion. As explained in detail in Appendix [Bj a negative 7 
at the phase transition in the thermodynamic limit im- 
plies that the inverse dielectric constant jumps to zero 
discontinuously. As for the inverse dielectric constant, 
we use a general fourth-order modulus to account for the 
two interactions, 

2 

T(ai,a,2) / , Ojfi^u^p^a'Jfiupcr-) (31) 

fI,V,p, <7—\ 



where 



lpu P a = (^) [(M,M,)(M p M a ) 



-3(M^ z M u , z M p , z M a , z }}. (32) 



The explicit derivation of Eqs. (31) and (32), is given in 
Appendix [B) 



B. Details of the Monte-Carlo simulations 

The Monte-Carlo updating scheme consists of trial 
moves for one or two particles at the same time, to 
a randomly chosen new location on the surface of the 
sphere. The change in the action Eq. (24) was calcu- 



lated and the move was accepted or rejected according to 
the Metropolis-Hastings algorithm 67 * 68 *. The trial moves 
were performed in three different ways. The first way 
was to move a single particle to a new random loca- 
tion uniformly over the total surface. The second way 
was to move a single particle to a new random loca- 
tion uniformly within some short distance, adjusted to 
yield a high acceptance rate. The last trial move was to 
move a nearest- neighbor pair of one z-particle and one 
^-particle together, to a random new location uniformly 
within some short distance, adjusted to yield a high ac- 
ceptance rate, and with a random new orientation. In 
order to straightforwardly ensure detailed balance, we 
additionally required the two particles to mutually be 
nearest-neighbors both in the old and the new config- 
uration. And to ensure ergodicity, the pair-move must 
be mixed with a number of single-particle moves. All of 
these moves were found to be essential in order to have 
fast thermalization as well as short autocorrelation times 
for the cases considered here. Pseudorandom numbers 
were generated by the Mersenne-Twister algorithm!^ and 
the sampled data were postprocessed using Ferrenberg- 
Swendsen reweighting technique J 70 * 71 1 



C. Results 

Motivated by its relevance to the fractional quantum 
Hall effect (in particular, the v = 1/2 MR state), we fo- 
cus on analyzing the screening properties of this system 
at Qi = 2 (M = 2pJ We also perform simulations in 
the neutral two-component Coulomb gas case at Qi =0 
(M = 0) in order to provide a check on the numerics, as 
well as for comparison with the Qi = 2 case. Further- 
more, the system is also studied for a number of values of 
the packing fraction, 77 to extract the screening properties 
in the low-density limit. 

For the two cases of Qi and the values of Q2 stud- 
ied below, the quantities e^ 1 and were found to be 
zero, within statistical uncertainty, and except for a small 
finite-size effect when system size TV, was small. Thus, we 
focus on the results for e^ 1 as this was the only term in 
Eq. (27) that contributed to the general inverse dielectric 
This means that screening properties 



constant, e 



-1 

(oi,a 2 ) ' 



of particles that interact with charges of both types, are 
determined by the charges of type 2, only. Note also that 
when = 0, the unconventional Coulomb plasma will 
screen test particles with charge of type 1, only. 

In Fig. [2J we plot e^ 1 in the relevant range of Q\ when 
the two-component neutral Coulomb gas (Qi = 0) is 
known to have a BKT transition. At small values of 
the system is in the screening phase where e^ 1 ~ 0. The 
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Figure 2: (Color online) Plot of the inverse dielectric constant 
for the model in Eq. |7| with Qi = and 1 < Q| < 
4.8. Results are presented for three different values of packing 
fraction 77 and three different values of system size N. 



Figure 3: (Color online) Plot of the inverse dielectric constant 
€22 for the model in Eq. ^ with Qi = 2 and 1 < Q\ < 
4.8. Results are presented for three different values of packing 
fraction 77 and three different values of system size N. 



reason for the ~ sign rather than an equal sign is that 
there is a mainly size-dependent offset from e^ 1 = 0, 
because perfect screening is not possible with a small 
number of charges. For large Q\ there is a phase in 
which charges of different components form tightly bound 
dipoles and the Coulomb gas turns into an insulator 
where e^ 1 ~ 1- Here, there is a mainly density-dependent 
offset from e^ 1 = 1 because the polarizability of the sys- 
tem increases with density, since the hard-core diameter 
d yields a minimum distance between the charges in the 
dipoles. The plot in Fig. [2] indeed shows that the charge- 
unbinding transition is dependent on the number of par- 
ticles in the system, as well as the size of the hard disk 
charges. When N increases, the onset of a finite value in 
moves to higher values of Q\. However, when we re- 
duce 77, the value of Q\ at onset of e^ 1 becomes smaller. 
Thus, this figure illustrates that understanding the be- 
havior in both limits N — >• 00 as well as 77 — » is not 
straightforward. 

In Fig. [3j results for the same case as in Fig. [2] are 
presented, but with Q 1 = 2. The results for Qi = and 
Qi = 2 are very similar, both qualitatively and quanti- 
tatively. Thus, the screening properties with respect to 
charge of type 2 of the unconventional Coulomb plasma 
when Qi = 2 are very similar to the well-studied two- 
component neutral Coulomb gas. 

To get a qualitative picture of the type 2 charge bind- 
ing of the unconventional plasma, three snapshots of the 
charge configuration when Qi = 2, 77 = 5 • 10 -4 , and 
N = 200 is given in Fig. [i] When Ql = 1, deep into the 
screening phase of the system (see Fig. [3|, most charges 
are free and only a small fraction of the charges may be 
said to form closely bound dipoles. At Q\ = 3, which 
is the relevant value for the Ising-type quantum Hall 
states, the system is closer to the unbinding transition 
and a larger fraction (though not all) of the particles are 



bound in dipoles. At Q\ = 5, deep in the type 2 insulat- 
ing region, all particles form closely bound dipoles and 
the ability to screen type 2 test charges is lost. 

Although it is clear from Figs. 2|3 that there is a tran- 
sition between a screening phase and an insulating phase, 
it is not easy to spot the transition point in the curves 
in these figures, which look rather smooth. Therefore, 
we must make some assumptions about the nature of the 
transition in order to identify it. 

For the case Qi =0, where the transition is known 
to be a BKT transition, it is natural to follow a method 
that was proposed in Ref. 72. At the BKT critical point, 
€22 scales logarithmically with N for large N. It taked 
the following finite-size scaling form: 



2 1 (iV)=e 2 - 2 1 (oo) 



ln(AT) + C_ 



(33) 



where e^ 2 1 (oo) is the value of ^(TV) when N — ^ 00 and C 
is an undetermined constant. Least-squares curve-fitting 
to Eq. (33) may be performed for various sizes N with 
C and €32 (00) as free parameters at fixed values of Q\. 
The critical point is then estimated as the value of Q\ 
which exhibits the best fit to Eq. (33). Additionally, for a 
BKT-transition, the value of e 2 " 2 1 ( 00) obtained at the best 
fit, must correspond with the universal jump condition, 



(52c e 22 1 (°°) = 4, cf. Eq. (11). Details of this procedure 



are given in Appendix [C] 

For Qi = 2, motivated by the similarity between 
Figs. |2j [3j we assume that the transition is also a BKT 
transition. We again look for the Q\ value at which the 
system best fits Eq. (33). Since we are able to find a 



value at which there is a very good fit to this form, we 
conclude that our assumption was justified. 

In Fig. [5j we present results for the criti- 
cal coupling Q\ c for four different densities 77 = 
0.0002,0.0005,0.001,0.002, for Q 1 = and Qi = 2. The 
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Figure 4: (Color online) Snapshots of the charge configuration 
at Q\ = 1, 3, 5 when Q 1 = 2, r] = 5 • 10" 4 and iV = 200. Red 
markers represent u>- particles, while blue markers represent 
^-particles. The marker diameters are about 5 times larger 
than hard disk diameter d. 




0.002 



Figure 5: (Color online) The critical value of Q\ found by 
curve fitting to Eq. (33) with two free parameters. Results 



are presented for four values of the packing fraction r\ and 
for Qi = and Q\ = 2. Fourteen system sizes in the range 
20 < N < 2000 have been used. 



results for Qi = reproduce the main features of the 
two-component Coulomb gas, namely that Q\ c = 4 when 
density is low and that Q\ c increases when density in- 
creases. These results also correspond well with earlier 
results in Refs. [Ml and [T5l When Q\ = 2, we find that 
the behavior of the critical temperature is very similar to 
the Qi = case, within statistical uncertainty. In addi- 
tion, in Fig. [6| results for the corresponding value of the 
parameter e^ 2 (oo) at the critical point is presented. The 
values for both Qi = and Qi = 2 are close to the uni- 
versal value of Q2c e 22 1 (°°) = ^ for the BKT-transition. 
Since the results for Q 1 = (the standard Coulomb- 
plasma BKT-transition case) and Qi = 2 are essentially 
the same, we suggest that the charge-unbinding transi- 
tion for the unconventional Coulomb plasma indeed is a 
BKT-transition in the sense that the type 2 inverse di- 
electric constant e^ 1 exhibits logarithmic finite-size scal- 
ing and a discontinuous jump with a universal value, as 
predicted by the BKT renormalization equations. 

As an additional verification of the discontinuous jump 
in the BKT-transition, we also study the fourth-order 
modulus 7( ai , a2 ), presented in Eqs. (31) and (32). As for 



the general inverse dielectric constant, we found that the 
only contributing term in the sum of Eq. (31 ) is the term 



with all indices equal to 2, 72222- Illustrating the typical 
behavior of this quantity, results for 72222 for a number 
of sizes when r] = 5 • 10 -4 and Qi = 2 are presented in 
Fig. [7J Typically, 72222 exhibits a dip at a value of the 
coupling that can be associated with the transition. As 
explained in Appendix [Bj a negative and finite dip in the 
limit when N — >> 00 signals the discontinuous jump in e^ 1 
that is a characteristic feature of a BKT-transition. To 
this end, the size of the dip in 72222 is plotted as a func- 
tion of inverse system size TV -1 in Fig. [8] in the case when 
77 = 5-10 -4 . The size of the dip |72222,min| decreases when 
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Figure 6: (Color online ) Th e universal jump value determined 
by curve fitting to Eq. ( |33| ) with two free parameters. Results 
are presented for four values of the packing fraction 77 and 
for Qi — and Qi = 2. Fourteen system sizes in the range 
20 < N < 2000 have been used. 
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Figure 8: (Color online) The size of the dip in the fourth- 
order modulus |72222,min| as a function of inverse system size 
A/" -1 . The packing fraction is r\ = 5 • 10 -4 , and results for 
Qi = and Qi = 2 are shown. The inset shows the results 
on a log-log scale. System sizes in the range 60 < N < 10000 
are used. 
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Figure 7: (Color online) The fourth-order modulus 72222 as 
a function of coupling Q\ for five different system sizes A/", 
when Qi = 2 and 77 = 5 • 10 -4 . 



AT increases towards the thermodynamic limit. However, 
assuming power-law dependence of |72222,min|, the pos- 
itive curvature in the log-log plot indicates a nonzero 
value of 1 72222, mini when N —> 00, verifying a discon- 
tinuous jump in , as expected for a BKT-transition. 
Again, we find that the results for Qi = 2 are very sim- 
ilar to Qi = 0. We also associate the coupling value of 
the minimum in the dip in 72222 with the critical point 
and the results are shown in Fig. [9] in the case when 
77 = 5 • 10 -4 . Clearly, the position of the dip moves to- 
wards higher values of Q\ when the system size increases. 
However, the evolution towards N~ x = is too slow to 
make a sharp determination of Q\ in this limit as also 
noted befor d 65 l 66 i With this method, we are not able to 
verify that Q\ Q ~ 4.4, as was found above in Fig. [5] for 
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Figure 9: (Color online) The coupling value at the minimum 
of the dip in the fourth-order modulus as a function of inverse 
system size A/" -1 . The packing fraction is r\ — 5 • 10" 4 , and 
results for Q\ — and Q\ — 2 are shown. The inset shows 
the results on a log-log scale. System sizes in the size 60 < 
N < 10000 are used. 



this density. 

By assuming a universal value of the discontinuous 
jump for a BKT-transition, we may determine the critical 



point of the BKT-transition using Eq. ( 33 ) with only one 
free parameter as described in Appendix C] The results 
The critical values of Q\ are very 
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are given in Fig. 

similar to what was obtained in Fig. |5j but are deter- 
mined with greater accuracy. For both cases, the critical 
point appears at higher Q\ when density increases. How- 
ever, Q\ c is systematically lower at Qi = 2 compared to 
Qi = 0. ' 

For the range of small densities that we have inves- 
tigated, the Monte-Carlo results for the unconventional 
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Figure 10: (Color online) The critical value of Q\ found by 
curve-fitting to Eq. ( |33[ ) with one free parameter. Results are 
presented for five values of the packing fraction r\ and for two 
values of Q\. 



Coulomb plasma with Qi =2 are rather conclusive. 
This plasma undergoes a charge-unbinding transition 
that should be regarded as a BKT-transition in the sense 
that the inverse dielectric constant of type 2 exhibits the 
well-established signatures of a BKT-transition. Specifi- 
cally, there is a density-dependent critical point Q\ c that 
separates a phase where particles of different species form 
bound pairs at high values of Q\ from a phase where par- 
ticles of different species are free at low values of Q\. For 
test particles carrying type 2 charge, the high-Q^ phase 
is unscreened, whereas the I0W-Q2 phase is screened. 

The results presented so far show that the behavior 
when Qi = and Qi = 2 are quite similar. However, 
in the phase with bounded dipoles, when charges of type 
2 are not screened, the cases Qi = and Q\ = 2 be- 
have rather differently. We first consider the case when 
Qi = 0. When charges are bound, this system consists 
of N dipoles that interact by dipole-dipole interactions. 
Consequently, these dipoles tend to form clusters with in- 
creased dipole strength, i.e. higher values of the coupling 



or the densitjff^^. In Fig. [Ill a snapshot of a Qi =0 
configuration with TV = 200, Q\ = 7 and rj = 2 • 10~ 3 
is shown, where some dipoles are seen to form clusters. 
In the case when Qi =2, the type 2 interactions are 
effectively reduced to dipole-dipole interactions, similar 
to the Qi = case. However, the logarithmic interac- 
tions of type 1 charges remain. Neglecting the weaker 
dipole-dipole interactions among dipoles of type two, the 
dipoles now essentially form elementary constituents with 
charge Qi interacting logarithmically. Effectively, the 
two-component unconventional plasma is reduced to a 
one-component plasma where the particles carry charge 
of type 1 and a (neutral) dipole of type 2. When Qi = 2 
this plasma is in the liquid state, i.e. the tightly bound 
dipoles do not form an ordered state with a broken trans- 




Figure 11: (Color online) Snapshots of the charge configura- 
tion at Qi = and Qi = 2 when Q\ = 7, 77 = 2 • 10" 3 , and 
N = 200. Red markers are u>-particles and blue markers are 
z-particles. The marker diameters are about 2.5 times larger 
than hard disk diameter d. 



interaction of type 1 charge will prevent the dipoles from 
forming clusters. A snapshot of the state with bounded 
dipoles when Qi = 2 is shown in Fig. 11 and the qual- 
itative difference from the case when Qi = is clearly 
seen. Quantitatively, this is seen by the behavior of e^~ 2 , 



presented in Fig. [12] When Qi =0, dipole-dipole in- 
teractions at short distances will reduce the fluctuations 
in the dipole moment resulting in a weakly increasing 
6*22 inside the bounded phase. On the other hand, when 
Qi = 2 the logarithmic interaction of type 1 charge will 
keep the dipoles at some distance from each other, thus 
the fluctuations of a dipole are not much restricted by 
the surrounding dipoles. Moreover, the strength of the 
dipoles increases with Q\ and a reduction in e^ 2 follows. 



lational or orientational symmetry. Also, the logarithmic The qualitative difference between the cases Q\ = and 
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Figure 12: (Color online) Plot of the type 2 inverse dielectric 
constant for Qi = and Qi = 2 with iV 100, 77 = 5 • 10~ 3 
in the range 3 < Q\ < 12. 

Qi = 2 is an effect due to the minimum separation of 
charges at finite density originating with the hard cores, 
and it will vanish in the limit 77 — ^ 0. 



IV. CONCLUSIONS 

We have shown that the unconventional Coulomb 
plasma analyzed in this paper, where particles can carry 
two distinct types of Coulombic charge, will screen test 
particles with charges of both types for the case most 
relevant for the plasma analogy of Ising-type fractional 
quantum Hall states, i.e. when there is one species of 
particles that carry type 1 charge Qi = 2 (M = 2) and 
type 2 charge Q2 = y/S and another species of particles 
that carry only type 2 charge — Qi- For test particles 
carrying both types of charge, screening will cease to oc- 
cur at Q\ = Q\ c « 4 in the limit of small density, when 
Qi = 2. For higher values of the system will continue 
to screen test particles that carry only type 1 charge, but 
will not be able to screen test particles with type 2 charge. 

One striking feature of these results is that Q\ c and the 
critical behavior at this point hardly seem to depend on 
Qi when density is small. This implies that the role of the 
type 1 interaction (which corresponds, in quantum Hall 
wavefunction language, to the Laughlin-Jastrow factor 
which accounts for the filling fraction of the system) is 
simply to maintain the ^-particles in a liquid state. Since 
its critical point is very far away, the type 1 interaction 
leads to a weak, smooth dependence on Q\. The physics 
in the transition at Q\ c is then dominated by the type 2 
interaction. We therefore conjecture that our results hold 
for all reasonable values of M - not only M = and 2, 
the cases which we have studied here, but also M = 1 
(which may be relevant to ultra-cold trapped bosons) and 
larger values of M, possibly all the way up to or near the 
critical value M c w 70, below which the one-component 



plasma of Eq. (p3l) is in the metallic phase^MD 
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Appendix A: Generalizing the inverse dielectric 
constant for multiple interactions 

In the unconventional plasma with two components 
that interact with two different Coulomb-like interac- 
tions, we are free to insert test particles that may interact 
with different charge strength through both interactions 
simultaneously. Here, we will generalize the inverse di- 
electric constant for such test particles. For consistency, 
we will also perform the derivation on the surface of a 
sphere by expanding in spherical harmonics. For a sim- 
ilar derivation, but with one interaction only and on a 
planar geometry, see Refs. I31|321 

When an external test charge field is inserted in the 
system, the free energy in the system will change accord- 
ing to the effective interaction among the test charges, 

AF[Sq] = [<m f dn'J2^^)U^(r-r')5q v (e',cj>'). 

(Al) 

Here, the effective interaction between charges of type \i 
and z/, is assumed to be of the form Uff = Uff(r • r'), 
Sq^O, <j)) is the test charge field for charges of type /i, and 
the integrations are over the solid angle dQ. To correctly 
model the test particles as carrying charge of different 
types, we write 

Sq fl (0^) = a fI Sqp(0^), (A2) 

where a M is a type dependent factor that accounts for 
the relative strength of charges of different types. For in- 
stance, the choice (ai,a 2 ) = (Q1/MQ2A) = (V 2 /3 M , 1) 
describes the test charges corresponding to quasiholes in 
the MR state, as given in Eq. (125) in Ref. 2 , which map 
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to particles in the plasma that carry charge Q\/2M = 
1/V2M of type 1 and charge Q 2 /2 = V^/2 of type 2. 
Moreover, in Eq. (A2) Sq is a common charge factor for 



all types such that a M Sq is the total charge of type \i 
carried by a test particle (which means that Sq = VS/2 
in the example above), and p(6, <j>) is the density field of 
the test particles. 

It is now convenient to expand the interaction and den- 
sity fields in spherical harmonics. The test particle den- 
sity field is expanded by 



p(M) = E E pTYrw,' 



(A3) 



1=0 m=-l 



where 

Yr(e,<j>) 



' {21 + - m)\ 
47r(Z + m)\ 



PF (cos 6)d 



imcf) 



(A4) 



and PJ n {x) are the associated Legendre polynomials. The 
coefficients are given by 



Pi 



(A5) 



The effective interaction is expanded by using the addi- 
tion theorem for spherical harmonics, 



uf^r •*') = £ airi S y r v> two 1 ' & F m 

1=0 m=-l 



(A6) 

Here, t are the Legendre coefficients of the interac- 
tion, given by 



rreff 
u fiu,l — 



2Z + 1 



/ d<9 sin flt/^(cose)fi (cos«), (A7) 
Jo 



where Pi(x) is the Legendre polynomial of order /. Now 
Eq. ( Al ) is written 



AF[5 q } = S q 2 Y J ^ l T, a ^ U f^ E 



1=0 



(A8) 

Hence, in the limit when the test charge field is infinites- 
imal, Sq —> 0, we find that 



d 2 F[Sq] 



dSq 2 



5q=0 1=0 



£^£v^£/r/>r. 



m— — l 



(A9) 

This derivative can also be calculated by inspection of 
the partition function of the system perturbed with the 
external test charge field. With F[Sq] = — In Z[Sq] and a 
potential energy on the form V[Sq] = Vo + V\[Sq) where 
Vo is the potential energy of the unperturbed system and 
V\ [Sq] is the contribution due to the test charge field, we 



find that 
d 2 Fl 



5q=0 




(MO) 



Here, we have also used that dF[Sq]/dSq\§ q =o = 0, and 
the brackets denote statistical average with respect to 
the unperturbed system. The test charges Sq^(6,(j)) will 
interact with each other as well as with the charge field 
g M (#,(/>). As for the test ch arge field, the charge field is 
expanded according to Eq. ( A3 ) to yield 



KlSq] = /dO f dO'£MM)+<MM)] 
J J fl 

xU(r-r')5q^8',<p') 

DO . I 

=E 2 m c/ <E^ E *«pr «i+V9pD- 

fi m— — l 

(AH) 



1=0 



whe re U (r • r) is the bare interaction, expanded by 
Eq. (A6) with coefficients U\. Performing the derivatives 
in Eq. (|A10|) yields 



dSq 2 



= E am Ul ^ a ^ vav ^ prpr 

Sq=0 1=0 W m=-l 

oo ^ OO ^ 

-EaTTT^EaTm^'E^ 

/=0 Z'=0 M,^ 

x E E 

m= — l m' = — l' 



We introduce the dielectric function e^j by 



(A12) 



(A13) 



and by comparing Eqs. (A9) and (A12), the inverse di- 
electric function is found to be 



\rn=-l J l'=0 

I V 

x E E pr^i?' (aw) 

m= — l m' = — l' 



Moreover, since the bare interaction is only dependent 
on the distance between the charges, U = U(r • r'), we 
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have that (q^q™/) = (<</)^<W, which yields 

/ / x -1 



— 1 e 



2tt 



-u. 



21 + 1 

= — L / 
I 

E^rvr (Ai5) 



Notice that even though there is no bare interaction be- 
tween charges of different type, there may be nonzero 
cross terms in Eq. (Al), as charges of different type are 



constrained to be together within the same particle. 



Appendix B: Fourth-order free energy derivative 



Additionally, the property that the bare interaction is 
distance dependent, only, yields an interaction U\ that is 
independent of m. Hence, the correlator (q™iq™*) must 
be m independent as well, (q™iq™*) = (q^iQui)- The 
dielectric function thus reads 



2tt 
2/ + 1 



Ul(q°M. 



(A16) 



The dielectric constant e^ v is now found in the long 
wavelength limit of the dielectric function. On a spherical 
surface, this corresponds to setting I = 1 in the dielectric 
function, i.e. e^ v = e^i. Thus, the dielectric constant 
is 



2tt 



(A17) 



So far, only a few assumptions are made regarding the 
bare interac tion U (r • r) and the charge field q^(6, </>). To 
apply Eq. (A17) for the system under consideration in 
this paper, we invo ke U (r • r) = — ln(l — r • r) to find 
Ui = 3/2 by Eq. (A7). Moreover, the charge field is 



modeled as point charges in a uniform background 



N 



5(6 - 0i)8{<t> - fa) 

i • n ; 



(A18) 



where q^ = — e /i,i)/(47r) is the uniform background 
ensuring charge neutrality for charges of type /i, e^^ is 
the charge of type \i in particle i and the sum is over all 
N p articles of the unperturbed system. Now, using Eq. 
( A5 ) , the actual coefficient of the charge field is found to 



be 



3 M n 



4tt R 



(A19) 



where M M = Y^iLi e /i,A is the total dipole moment for 
charges of type \i. Finally, by inserting these results in 
Eq. (A17), the inverse dielectric constant is found to be 



J \1V 



(A20) 



where (M^M^z) = (M M • M^)/3 by assuming isotropy. 

When there are test charges with multiple interactions, 
there are multiple contributions to the change in free en- 
ergy as seen in Eq. ( Al ). To account for all contributions 



to the increase in free energy, we construct a generalized 
dielectric constant by 



"(ai,a 2 ,...) 



€ jJLV a ^ 



(A21) 



In Ref.|65]a method of verifying the discontinuous char- 
acter of the BKT-transition was introduced, by examin- 
ing a higher-order term in the free energy expansion in 
the XY-model when the system is perturbed with an in- 
finitesimal phase twist. Similarly, in Ref. [66j the method 
was applied in a two-dimensional logarithmic plasma. 
Here, we show that the same idea also applies when we 
perturb a logarithmic Coulomb plasma on a spherical sur- 
face with an infinitesimal test charge field with multiple 
types of Coulomb interactions. 

Consider a system with particles interacting with dif- 
ferent charges of multiple types, as previously described. 
We now choose to perturb this system with a neutral dis- 
tribution of test charge of multiple types, which has the 
form Sq fI (6) = a^Sq cos(0), i.e. a similar test particle den- 
sity field as given in Eq. ( |A2[ ) but with = ^4tt/3 being 
the only nonzero coefficient in the spherical harmonics 
expansion. This is a convenient choice because it corre- 
sponds to the most long-waved nonuniform test charge 
configuration on the surface of a sphere, and hence, the 
prefactor of the second-order term in the free energy ex- 
pansion will be proportional to the inverse dielectric con- 
stant, as we will see below. 

The test charges yield a contribution to the po tentia l 
energy as given by the I = 1 and m = term in Eq. ( All ), 



Vi[Sq] = (<# ,i + 



(Bl) 



We now consider how the system responds to the test 
charges by a Taylor expansion of the free energy in the 
test charge field around 5q = 0, 



AF[5q] 



dF[Sq] 



Sq 



d 2 F[5q] 



5q=0 

d 3 F[Sq] 



Sq=0 



5£ 

2! 



d5q 3 



Sq=0 



8tf_ 
3! 



d A F[8q] 



d5q 4 



Sq=0 



Sq*_ 
4! 



(B2) 



The change in the free energy AF[Sq] must be invariant 



to Sq^O) 



-Sq^O), and hence, all odd-order deriva- 



tives in Eq. ( |B2| ) a re ze ro. From Appendix |A| (see Eqs. 
( |A12| ), ( |A17 ) and ( A21| )), the second-order free energy 
derivative is found to be 



d 2 F[5q] 



8tt 



5q=0 



(B3) 
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The fourth-order derivative is 
d A F\5q] n I [ dV^Sq] 



dSq 4 



5q=0 




Sq=0 j 



dV^Sq] 



Sq=0 j 



V 3 

x - ( B4 ) 

where brackets denote a statistical averag e wi th res pect 
to the unpe rturbed action. Inserting Eqs. (B3) and ( |B4| ) 
in Eq. (|B2| yields 



AF[5q] = f (p?) 2 /7x 



e (ai,a 2 ,...) 2! 
+ 7(ai,a 2 ,...) + ' 



where 



and 



T(ai,a 2 ,...) — ^ ^ ®>iL®>vQ>p®>(j^iiLvp(ji 



(B5) 



(B6) 



/47r TT \ 3 (p?) 2 



Now, insert ing = ^4tt/3 and assuming the charge 
field in Eq. (A18) and a logarithmic bare interaction, 
C7i = 3/2, yields 

7^ = (^) 2 [< M m M -> (M p M a ) 

- 3 (M^ z M v , z M p , z M a , z ) } , (B8) 

where {M^^ z M ViZ ) = (M M • M^)/3 by assuming isotropy. 

1. Stability argument 

When Sq = 0, the free energy of the system has a 
global minimum, and hence, the right-hand side of Eq. 
(B5) must be greater or equal to zero. Now, if 7( ai ,a 2 ,...) 
approaches a nonzero negative value at the critical point 
in the thermodynamical limit, the general inverse dielec- 
tric constant must simultaneously have a nonzero posi- 
tive value for the ground state to be stable. However, 



since e 



(ai,a 2 ,...) 



~(ai,a 2 ,...) 



in the screening phase, it follows that 
must exhibit a discontinuous jump at the criti- 
cal point. Hence, investigation of 7( ai ,a 2 ,...) ma y be used 
to verify a discontinuity in the inverse dielectric constant, 
which is a necessary requirement for observing a BKT- 
transition. 



Appendix C: The finite-size scaling relation 

The finite-size scaling relation of the BKT-transition 
has been used throughout this article to verify the uni- 
versal jump in e^2 and to provide estimates for the criti- 
cal coupling Q\ c . Here, some details to the curve fitting 
procedure and the goodness of fit measure are presented. 



1. Two free parameters 

Least-squares curve fitting of the Monte-Carlo results 
for to Eq. (33) may b e performe d with both e^~ 2 (oo) 
and C as free parameter a 28 l 3Q l 72 l 73 [ If the transition is 
of the BKT-type, a good fit to Eq. (|33|) should be ob- 



tained at the critical point. In addition, when (^(oo) is 
free, no a priori assumption on the value of the univer- 
sal jump is made, thus a resulting value of e^ 2 (oo) that 
corresponds to the universal jump of the BKT-transition 
should be obtained. However, with two free parameters, 
higher quality of the Monte-Carlo statistics is required 
to single out when they system is closely obeying the 
behavior of Eq. ( 33 ) . 



We have employed the Marquardt-Levenberg algo- 
rithm minimizing y 2 to the nonlinear fitting function in 
Eq. ( |33] ). Specifically, x 2 is the sum of squared weighted 
residuals, 



(CI) 



where n is the number of system sizes iVj, e 22 1 Ar . is the 

value of the inverse dielectric constant e^ 1 obtained from 
the Monte-Carlo simulation at system size N{, and 
is the corresponding error. For a good fit, we expect the 



weight-normalized residuals, Y{ = (e^ 1 



22, Ni ' 



c 22 



!(Ni))/a Ni 



to be Gaussian-distributed with mean n(Yi) = and 
variance cr 2 (F z ) = 1. Thus, to measure the goodness 
of the fit, we use the Anderson-Darling test statistic A 2 
for the data set Y\ to arise from a normal distribution 
with fi(Yi) = and a 2 (Yi) = 1: 

n 

2=1 

(C2) 

where &(Y) is the standard normal cumulative distribu- 
tion function and where the data set Yi is ordered from 
low to high values. A smaller value of A 2 essentially 
means a better fit between the data and the fit function. 

To illustrate the method, Monte-Carlo results for e 22 
at fourteen different system sizes and the corresponding 
curve-fit according to Eq. (33) are given in Fig 
three different values of Q% 
0. Clearly, at Q 



Here, rj = 2 • 10" 
f\ = u. ^leciriy, ciu Q>\ = 4.933, the fit between 
data and the fit function is better than for the two other 



for 
and 
the 



cases. Moreover, in Fig. 14 the corresponding results for 
the goodness of fit parameter as well as the results for 
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Figure 13: (Color online) Plot of the size-dependence in the 
inverse dielectric constant e^ 2 1 (A/') for fourteen different sys- 
tem sizes in the range 20 < N < 2000 at three different values 
of the coupling Q\. The best fit according to the fit function 
in Eq. (33) with two free parameters, is given as the corre- 



sponding solid line in all three cases. The packing fraction is 
77 = 2 • 10" 3 and Qi 0. 



the parameter e^ 2 1 (oo) as a function of Q 2 are shown. 
Indeed, the minimum in A 2 indicates a critical region 
where the data seem to follow the logarithmic finite size 
scaling of e^ 1 given in Eq. (33). Also note that this region 
coincides with a value of yp^H 00 ) c l° se to the universal 
jump value of 4. With the minimum of A 2 as a measure 
of the critical point and with error estimates obtained by 
the Jackknife method, we find that Q 2 2 c = 4.933 ± 0.012 

and that Ql^ii 00 ) = 3 - 941 ± °- 023 > less than 2 % off 
the universal number. The results in Figs. [5] and [6] are 
found by repeating this procedure for different values of 
77 and Qi. 



2. One free parameter 

The procedure described in detail above with two free 
parameters, may be performed with a fixed value of 
e^ 2 1 (oo) = 4Q2 c an d with C as the only free parame- 
ter. If the transition is of the BKT-type, a good fit to 
Eq. ( [33] ) should be obtained at the critical point. This is 
a rather well-used m ethod to d etermine the critical point 
of a BKT-transitio ri 29 ' 72 ' 74 ' 75 4 With only one free pa- 
rameter, Q\ c will be determined with greater accuracy 



compared to the case when there are two free parameters. 



3. Remarks 

Refs. [28] and [30] used x 2 as a goodness of fit parame- 
ter. We also tried this, and the results for the critical cou- 
pling as well as the corresponding parameter e^ 1 (00) were 

consistent with A 2 results within statistical uncertainty. 

50 1 , , , 1 4.2 




Figure 14: (Color online) Plot of the goodness of fit parameter 
A 2 and the corresponding free parameter e^" 1 (oo) obtained 
when curve fitting to the critical finite-size relation given in 
Eq. ( |33| ). The results are given as a function of Q\. System 
sizes N, and r\ and Qi are the same as in Fig. |13| Error 
estimates are obtained by the Jackknife method. 



However, we found that error estimates were clearly un- 
derestimated with x 2 , probably due to over- fitting. 

The parameter C in the finite-size scaling relation (Eq. 
([33]) ) is density dependent!^. Specifically, C increases 
when rj decreases. Hence, at the critical point, the finite- 
size scaling slows down when rj is lowered. Therefore, 
larger systems N or better statistics are required to re- 
solve the critical scaling when 77 is small. In partic- 
ular, curve fitting to Eq. (33) was also performed for 
10 -5 in addition to the densities presented in 



T] = 5 

Figs. [5] and [6] However, in this case the statistics were 
not good enough to resolve a clear minimum in A 2 . Also 
note that there are higher-order correction^ to Eq. (33) 
that are not taken into account in this work. 
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